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ABSTRACT 

Context. The deceleration mechanisms for relativistic jets in active galactic nuclei remain an open question, and in 
this paper we propose a model which could explain sudden jet deceleration, invoking density discontinuities. This is 
particularly motivated by recent indications from HYbrid MOrphology Radio Sources, suggesting that in some case 
Fanaroff-Riley classification is induced by variations in density of the external medium. 

Aims. Exploiting high resolution, numerical simulations, we demonstrate that for both high and low energy jets, always 
at high Lorentz factor, a transition to a higher density environment can cause a significant fraction of the directed jet 
energy to be lost on reflection. This can explain how one-sided jet deceleration and a transition to FR I type can occur 
in HYbrid MOrphology Radio Sources, which start as FR II (and remain so on the other side). 

Methods. For that purpose, we implemented in the relativistic hydrodynamic grid-adaptive AMRVAC code, the Synge- 
type equation of state introduced in the general polytropic case by Meliani et al. (2004). To demonstrate its accuracy, 
we set up various test problems in appendix, which we compare to exact solutions that we calculate as well. We use the 
code to analyse the deceleration of jets in FR II/FR I radio galaxies, following them at high resolution across several 
hundreds of jet beam radii. 

Results. We present results for 10 model computations, varying the inlet Lorentz factor from 10 to 20, including uniform 
or decreasing density profiles, and allowing for cylindrical versus conical jet models. As long as the jet propagates 
through uniform media, we find that the density contrast sets most of the propagation characteristics, fully consistent 
with previous modeling efi^orts. When the jet runs into a denser medium, we find a clear distinction in the decelaration 
of high energy jets depending on the encountered density jump. For fairly high density contrast, the jet becomes 
destabilised and compressed, decelerates strongly (up to subrelativistic speeds) and can form knots. If the density 
contrast is too weak, the high energy jets continue with FR II characteristics. The trend is similar for the low energy 
jet models, which start as underdense jets from the outset, and decelerate by entrainment in the lower region as well. 
We point out differences that are found between cylindrical and conical jet models, together with dynamical details 
like the Richtmyer-Meshkov instabilities developing at the original contact interface. 

Key words. ISM: jets and outfiows - Galaxies: jets - methods: numerical, relativity 



1. Introduction 

Accretion disks surrounding black holes, jets found in as- 
sociation with compact objects, and Gamma Ray Bursts 
(GRBs) all represent violent astrophysical phenomena. 
They are associated with relativistic flows, both with re- 
spect to the occuring velocities and to their prevailing 
equation of state. The jets from Active Galactic Nuclei 
(AGN) and in GRBs are accelerated in a short distance 
to reach a hig h Lorentz factor: typical values are 7 '--^ 5 — 30 
KeUermann et a l. 2004) for AGNs and 7 > 100 for GRBs 
Woods fc: Loeb||1 995 Sa ri fc Piran||1995 1. In this accelera- 



tion phase, situated at the base of the jet, it is believed that 
jet energy is dominated by thermal energy and Poynting 
flux, and that a fraction of these energies contributes to 
the efficient acceleration of the flow. Thereby, the relativis- 
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tic fluid changes its state from relativistic, corresponding 
to an effective polytropic index F = 4/3, to classical (poly- 
tropic index F = 5/3) when the thermal energy is con- 
verted to kinetic energy. Also further in the jet paths, where 
the jets interact with the surrounding medium, a fraction 
of the directed kinetic energy is converted to thermal en- 
ergy at the shock fronts formed. According to the prevailing 
Lorentz factor of the jet, these shocks could be relativistic 
and therefore very efficient to convert kinetic to thermal 
energy, or could be Newtonian and only weakly efficient in 
compression. To investigate the occuring relativistic ffows, 
a realistic equation of state must therefore be adopted, to 
handle both classical as well as relativistic temperature 
variations in space and time. For that purpose, we im- 
plemented in the relativi stic (magneto) hydrodynamic grid- 
adaptive AMRVAC code ([M eliani et al. 2007||Keppens et al. 
2003 van der Hoist fc Keppens^200f ), the Synge-type equa- 



2 
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tion of state introduced for a general polytropic case as in 
Meliani et al. (2004) and as applied in the adiabatic case by 
Mignone & McKinney (2007). To demonstrate its accuracy, 
we set up various stringent test problems in an appendix, 
which we compare to exact solutions for Riemann problems, 
that we calculate as well. The latter include Riemann prob- 
lems at Lorentz factors of order 7 « 100, which represent 
extreme values relevant for Gamma Ray Burst flows. 

The development of relativistic numerical hydrody- 
namic codes can help us understand the physics of as- 
trophysical jet propagation. In the last decade, significant 
progress was made in numerical special relativistic hydro- 
dynamic (HD) and magnetohydrodynamic codes. Various 
authors worked on the development of conservative shock- 
capturing schemes which use either exact or approximate 
Riemann solver based methods for relativist ic hydrodynam- 
Mellema 
Zanna & 



propagation in uniform overden se media and use an equa- 
tion of state for a pair-plasma. Scheck et al.' (2002| study 



ics ( Eulderink 



et al.,|1999 



Mellema| |1994| [Font et aT] p94| [A% 
B 



ucciantini 



20021 IMignone fc 



Bodo| [20051 (for a review see [Marti fc Muller| ^2003| )) 
i'he study of relativistic hydrodynamic fluids benefits also 
from using spatial and te mporal adaptive techniques, or 
Adaptive Mesh Refinement ([Dunca n & Hughes 1994; Zhang 
fc MacFadyen|2006[ [Meliani et al.[[2007( .Wang et al.|2007| 



The various numerical simulations usually involve a simpli- 
fied equation of state (EOS) with a constant polytropic in- 
dex. Notable exceptions with more realistic EOS treatments 



are found in Mignone et al. 


([20051); 


Mignone & McKinney 


(20071; 


Perucho &: Marti 1 


2007|). We present in the ap- 



pendix to this paper the required extension of the AMRVAC 



code (Meliani et al. 2007 Keppens et al. 2003 van der Hoist 
[fc Keppens 200 7 1 with the more re alistic polytropic EOS 
introduced by Meliani et al. (20041 which is based on the 
Syn ge ga s equati on, the relativis tically correct perfect gas 
law ( [Syn ge 1957; Mathews[[1971 ). The equations which we 
solve, and the schemes used, are also mentioned there. 

Many previous investigations of relativistic jet prop- 
agation through the interstellar medium (IS M) in radio 
galaxies concentrate on uniform ISM conditions (Duncan &I 
Hughes., 1994.; Marti et al. 1997, ; Komissarov fc Falle 1998[ 



Aloy et aL]T999[) . These investigations contributed to our 
understanding of the jet deceleration, where one then dis- 
tinguishes dynamics for Faranoff-Riley type I (FR I) and 
FR II radio galaxies, according to the power of the jet and 
hence the accretion rate in their galactic center. In the FR I 
radio galaxies, the associated jets are relativistic on parsec 
scale and sub-relativistic on kiloparsec-scales, so that jet 
deceleration, and thu s energy redistributio ns, must happen 
on kiloparsec scales (Hardcastle et al. 2005 1. Various studies 



have looked into possible deceleration processes with both 
non-re lativ istic and relativistic H D codes. IHooda fc Wiital 
([1996') and' Hooda fc Wiita[ ( [l998[ ) investigated with a clas- 
sical HD code the propagation of a 3D jet through an inter- 
face between dense ISM and less dense intercluster medium. 
They found that for such interface marking a density de- 
crease, the jet does not decelerate when crossing it, and that 
the deceleration mostly happens in the lower region (dense 
ISM). Moreover, in their study the jets are not defiected 
at the interface ISM /ICM. [Norman et"al] ([l988 ) also use a 
classical code to study the sudden deceleration of jets, when 
the jet crosses a shock wave in the external medium. Other 
works investigate the interaction of jets with dense clouds 
in both non -relativistic hydrodynamic si mulations ([Saxton 
et al. [[20051) and relativistic sin i ulatio ns ( Choi et al.[[2007p . 
Duncan, Hughes & Opperman ( 1996 1 study relativistic jet 



the influence of the matter composition of a high jet energy 
jet on its interaction with the external medium. They inves- 
tigate th e two extreme case s of pu re leptonic and baryonic 
plasmas. Perucho & Marti (20071 examined the propaga- 
tion of the low energy jet of the speciflc FR I radio source 3C 
31 using an elaborate hydrodynamics model, with an equa- 
tion of state that distinguishes the contribution of leptons 
and baryons to density and pressure. Most recently, [Rossi[ 
et al. ( 2008 1 investigate the propagation of jets in uniform 
overdense media in 3D, using a realistic synge EOS. 

In this paper, we investigate a scenario of relativistic jet 
propagation through an ISM with a sudden jump in den- 
sity. Since jets propagate over enormous distances, it is in- 
evitable that they encounter density jumps in the external 
medium. As we assume that matter in the inner part of the 
radio galaxy has been cleared away during the evolution of 
the radio galaxy, we typically follow jet dynamics through 
a lighter medium at flrst, which then suddenly changes far 
away to a denser medium. In a sequence of model runs, we 
vary the jet and medium properties to cover both high and 
low jet beam kinetic luminosities, straight as well as coni- 
cal jets (varying the jet opening angle), and allow for either 
uniform or decreasing density proflles. To put our work into 
perspective and to appreciate its relation to previous works 
better, we compiled in Table [T] the most relevant parame- 
ters and simulation specifics for selected relativistic hydro 
models. The most novel assumption is the transition from 
low to high density across a contact interface (across which 
the pressure is nec essarily constant ) , whic h is different from 
the work done by Hooda & Wiita ( 1998 ) where they use a 
classical code witn a density decreas e across the contact , 
and also different from the study by Loken et al. (1993), 



where they use a classical code to look into jet propagation 
through a pressure wall. We explore the influence of a sud- 
den density jump in ISM on jet propagation, stability, and 
formation of the bow shock. We aim to better understand 
the jet efficiency in transporting energy and mass from the 
central AGN regions to the denser ISM at large distances, 
where the jet cocoon gets formed. Note also from Table [T] 
that our models extend the previous works most notably 
by typically covering much larger distances than previously 
studied (as we need to obtain representative endstate mor- 
phologies in both lower and upper media), and are invari- 
ably in the high Lorentz factor regime, combined with a 
high resolution through the jet beam. The high energy jets 
we consider start out as denser than their immediate sur- 
roundings (a likely property deduced from all common jet 
launch scenarios, but previously ignored by all jet simula- 
tions to date). Moreover, the density contrasts we inves- 
tigate between jet and outer medium are much more rea- 
sonable than the 5 orders of magn itude density differences 
previously used by [Scheck et al. ((2002) . In what follows, 
we first motivate our model assumptions and the simula- 
tion setup in Section [2[ We discuss our main findings in 
Section [3| 



2. Relativistic jet propagation and deceleration 

2.1. Motivation 



The radio loud Fanaroff-Riley galaxies (Fanaroff & Riley 
1974) have extensively been studied in the last decade. 
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Table 1. The most re levant characteristics and param eters for various se lected relativistic hydrodyn amic s imulations in 



IW), 



the literatur e: DH94 (Duncan fc Hughes"1994') , M97 (|Martf et al.| |1997'), H98 CHarde e et^ 

p99 ) S02 (Scheck et al. 2002), PM07 ( Pcrucho fc MartilpOOTf T^S (Rossi ct al. 2008| . Our simmauons are maicaiea 
by MKG. We mention the type of the simulation 2D or 3D and use of AMR, jet beam kinetic luminosity Lb (when 
available), 7b Lorentz factor of the beam, Mb relativistic Mach number, rj density ratio ("d" means density variation), 
9 open angle of the jet, i?b/''ccn grid cells through jet beam, type of EOS ("Lep" means Leptonic, "Bar" Baryonic), 
endtime of the simulation in units of light crossing time of the jet beam radius. 



Paper 


case 


Lb 


7b 


Mb 


pmcdium/ Pjct 


Ob 


-Rb/fcoll 


EOS 


Size in _Rb (and time) 


DH94/R99 


A (2D, AMR) 




1.048 


6 


10.0 





24 


5/3 


10 X 41.67 (-) 


DH94/R99 


B (2D, AMR) 




5.0 


8 


10.0 





24 


5/3 


16.67 X 41.67 (-) 


DH94/R99 


C (2D, AMR) 




10.0 


8 


10.0 





24 


5/3 


16.67 X 41.67 (-) 


DH94/R99 


D (2D, AMR) 




10.0 


15 


10.0 





24 


4/3 


16.67 X 41.67 (-) 


M97 


Al (2D) 




7.1 


9.97 


10^ 





20 


4/3 


10.5 X 50 (50) 


M97 


A2 (2D) 




22.37 


31.86 


10^ 





20 


4/3 


10.5 X 50 (50) 


M97 


al (2D) 




7.1 


9.97 


1.0 





20 


4/3 


10.5 X 25 (25) 


M97 


a2 (2D) 




7.1 


9.97 


10.0 





20 


4/3 


10.5 X 25 (25) 


M97 


Bl (2D) 




7.1 


41.95 


10^ 





20 


4/3 


10.5 X 50 (50) 


M97 


B2 (2D) 




22.37 


132.32 


10^ 





20 


4/3 


10.5 X 50 (50) 


M97 


bl (2D) 


_ 


2.29 


13.61 


10.0 





20 


4/3 


10.5 X 25 (25) 


M97 


b2 (2D) 


_ 


7.1 


41.95 


10.0 





20 


4/3 


10.5 X 25 (25) 


M97 


CI (2D) 


_ 


2.29 


13.61 


10^ 





20 


5/3 


10.5 X 50 (50) 


M97 


C2 (2D) 


_ 


7.1 


41.95 


10^ 





20 


5/3 


10.5 X 50 (50) 


M97 


C3 (2D) 


_ 


22.37 


132.32 


10^ 





20 


5/3 


10.5 X 50 (50) 


M97 


cl (2D) 


- 


2.29 


13.61 


10.0 





20 


5/3 


10.5 X 25 (25) 


M97 


c2 (2D) 


- 


7.1 


41.95 


10.0 





20 


5/3 


10.5 X 25 (25) 


H98 


B (2D, AMR) 


- 


5.52 


8.87 


10.0 





24 


5/3 


16.67 X 41.67 (140.0) 


H98 


C (2D, AMR) 


- 


14.35 


11.52 


10.0 





24 


5/3 


16.67 X 41.67 (21.43) 


H98 


D (2D, AMR) 


- 


10.0 


15.0 


10.0 





24 


4/3 


16.67 X 41.67 (35.273) 


H98 


E (2D, AMR) 


- 


2.55 


8.53 


10.0 





24 


5/3 


16.67 X 41.67 (140.0) 


S02 


A (2D) 




6.62 


9.24 


lO'' 





6 


Lep 


200 X 500 (5950) 


S02 


B (2D) 




6.62 


14.33 


lO"" 





6 


Lep 


200 X 500 (5400) 


S02 


C (2D) 


10"*' 


7.95 


130.14 


lO"" 





6 


Bar 


200 X 500 (5300) 


PM07 


A (2D) 




2.0 


4.75 







16 


Lep/Bar 


200 X 450 (37231) 


R08 


A (3D) 


- 


10.0 


28.3 


10^ 





20 


Synge 


50 X 150 X 50 (240) 


R08 


B (3D) 




10.0 


28.3 


10 





20 


Synge 


60 X 75 X 60 (760) 


R08 


C (3D) 




10.0 


28.3 


10* 





12 


Synge 


50 X 75 X 50 (760) 


R08 


D (3D) 




10.0 


300 


10* 





20 


Synge 


50 X 150 X 50 (760) 


R08 


E (3D) 




10.0 


300 


10^ 





12 


Synge 


24 X 200 X 24 (150) 


MKG 


A (2D, AMR) 


10*« 


20.0 


1200 


0.1496-4.687 





64 


Synge 


10 X 400 (380) 


MKG 


B (2D, AMR) 


^q46 


20.0 


1200 


0.1496-671.22 





128 


Synge 


10 X 400 (900) 


MKG 


C (2D, AMR) 


10« 


10.0 


39 


36.52-203.66 "d" 





76 


Synge 


10 X 400 (820) 


MKG 


D (2D, AMR) 


10« 


10.0 


39 


36.52-148.15 "d" 


1 


76 


Synge 


10 X 400 (820) 


MKG 


E (2D, AMR) 


10« 


20.0 


35 


146-847.46 "d" 





76 


Synge 


40 X 400 (820) 


MKG 


F (2D, AMR) 


10« 


20.0 


35 


146-645.16 "d" 


1 


76 


Synge 


40 X 400 (820) 


MKG 


G (2D, AMR) 


^q46 


10.0 


1300 


0.033-0.495 "d" 





144 


Synge 


40 X 400 (380) 


MKG 


H (2D, AMR) 


10« 


10.0 


1300 


0.033-30.6748 "d" 


1 


76 


Synge 


40 X 400 (800) 


MKG 


I (2D, AMR) 


10"« 


10.0 


1300 


0.033-0.306748 "d" 


1 


76 


Synge 


40 X 400 (480) 


MKG 


J (2D, AMR) 


^q46 


20.0 


1200 


0.1496-12.95 "d" 


1 


76 


Synge 


40 X 400 (480) 



because their jets show intriguing behaviour. In fact, ra- 
dio loud galaxies are grouped in two main classes accord- 
ing to their radio map morphology. A typical FR I is 
brighter near the center and fades out towards the edge, 
whereas FR II a re brightest at the edge s and fainter to- 
ward the center. Fanaroff fc Riley ( 1974[ ) discovered that 
the break in FR I/FR II occurs around the radio luminos- 
ity of PiTsMHz ~ 10^^WHz~^sr~^, with almost all sources 
below the break value being of type FR I. In contrast to FR 
II galaxies where the jet re mains relativistic and n arrow in 
all scales, the jet in FR I ( Giovannini et al. 2005) is rela- 
tivistic at the parsec-scale ( Bridle|1992 ) and in many cases 
subrelativistic and diffuse in kparsec-scales (|Giovannini et 



the decrease of Doppler b eaming in an intr insically sym- 



metrical, decelerating jet (Laing et al. 19991. However, in 
some FR I the structure of the jet in the kiloparsec scale 
appears more complicated, with an inner spine which re- 
mains relativistic and a n outer shell that d ecelerates and 
becomes sub-relativistic ( Canvin et al.|[2005 ). 

It is commonly accepted that the morphological differ- 
ences between FR I/FR II arise because of differences in the 
physical conditions of jet interaction with its environment. 
However, it is still under debate whether the observed de- 
celeration in FR I jets is related to the mechanism of jet 
launching and thus to the properties of the central engin e 



in AGN (GhiseUini fc Celotti 



2001 Kaiser fc Best 2007), 



al. 2001 1 . The FR I jet shows side to side asymmetries which or rather related to the properties of the external mediuni 



decrease with distance from the central engine, caused by the host galaxy (|Zirbel|[l997| and the circumgalactic gas 



4 
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jPe" Young 



19931 iKaiser & Alexander||1997|. Some authors 

wn 



g||i99a . 

Id et al'l ( |2007| ) claim that the b"^ l/FR II dichotomy 
is caused by the combination of a central engine factor, 
with a contribution of the external medium. Observational 
evidence for the external medium influence comes from 
the existence of a group of peculiar "HYbr id MOrphology 
Radio Source" (HYMORS), pointed out by iGobal-Ktishna' 
I&: Witta (2002j . These radio sources appear to have an FR 
11 type on one side and an FR I type diffuse radio lobe 
on the other side of the active nucleus. This supports the 
idea that the different Fanaroff- Riley morphologies are at- 
tributed in some cases to the properties of the ambient me- 
dia, since for HYMORS, similar jets (power, composition, 
Lorentz f actor) likely emerge from the central engine on 
each side. Heywood, Blundell & Rawhngs ( 2007[ ) analysed 
radio images for a set of quasars 7C (radiosource galaxies) 
and found that some FR I seem to have a radio luminosity 
exceeding the original FR I/FR II dividing line. They also 
confirmed the existence of HYMORS and suggested that 
these sources have high power jets (typical for FR II) and 
they yield FR I or FR II morphologies according to their 
external medium. In a rarefied medium, the jet is 'lami- 
nar' and remains collimated all the way to the intergalactic 
medium (IGM) where the hotspot forms (the impact site 
of the jet in the IGM). This gives an FR II jet morphol- 
ogy with the typical "lobe" . In dense or clumpy medium, 
the jet interacts strongly with its external medium and dis- 
sipates its energy more gradually as it propagates in this 
dense medium. This gives the FR I jet morphology with 
the typical "plume" at large spatial scales. 

We point out several concrete examples next. The bright 
one-sided (FR I-like), diffuse jet in the radio galaxy 3C 321 
(which is an FR II radio source) is an excellent example of 
a HYMORS. Evans et al. (2008) argue that the one-sided 
diffuse jet is the result of an interaction of the jet with the 
companion galaxy. The fact that in these objects, both sides 
of the jets are relativistic and stable at the parsec scale, 
and the difference appears at the large scale, means that 
the variation in the external medium must occur at some 
distance from the central engine. This is what we assume in 
our models. The radio galaxy Gen A shows also a difference 
in the radio morphology between the two sides of the jet 
as it propagates on kiloparsec scales, with edge-brightened 
lobes (FR I-like) on one side, and on t he other side a cen - 
tral (fine structured) lobe (FR Il-like) (|Kraft et al.||2003l 



Yet another example is the powerful radio source Hercules 
A (3G 348) which has a jet kinetic luminosity lO^^ergs/s 



(McNa mara et al.||2005 l and exhibits a mixed FR I/FR 
II morphology ( Sadun fc Morrison^ 2002) . These observa- 
tions confirm the idea that in Gen A and Hercules A (3G 
348), the external medium plays a key role in the appear- 
ance of the jet and its dynamics. There is also evidence for 
disruptions or variations of radio morphologies induced by 
inhomogeneous medium. For example th e evident hierarchi - 
cal structure of M87 in its radio image ( Owen et al.||2000 1, 
with the possibility of the existence of two haios: an mner 
halo that could be with more porous s tructure, while th e 
jet interacts mainly with the outer halo ( Owen et al.|2000" l. 
In conclusion, som e FR I jets propagate through clumpy 
( Groft et al. 2006 1 or dense media, evidently encountering 
sudden density changes. This is known to give rise to a 
strong interaction, and the jet loses its energy by entrain- 
ment and di ffusion ( De Young|1993 1996 Rosen fc Hardee 
2000 [ 2006 1 and also forms knots along the jet ( Owen et 



al. 19891. We now detail our new, more elaborate model, 
with which we aim to explain observations of jets in radio 
galaxies that are relativistic on small scale and decelerate 
to sub-relativistic velocities going to the large scale. 

2.2. Model 

To make the problem more tractable, we assume axisym- 
metry and we neglect the influence of the magnetic field in 
the dynamics. In our scenario, we assume the existence of 
a density jump in the host galaxy or in circumgalactic gas 
which could be responsible for the jet deceleration and knot 
formation in the asymptotic region. Since the jets travel 
for enormous distances, it is inevitable that they encounter 
various sudden transitions of interstellar medium proper- 
ties. These could be traveling shock fronts, more gradu- 
ally varying background variations as one traverses regions 
of differing gravitational potential, or contact discontinu- 
ities indicative of boundaries between varying regions of 
influence. We concentrate on the latter, representing den- 
sity (and entropy) changes across which total force balance 
holds (uniform pressure), as these are invariably found in 
any kind of hydrodynamic interaction involving winds, out- 
flows, etc. 

We model jet propagation through two distinct media. 
In the first part, a low and uniform density is assumed, such 
that in this region the jet-external medium interaction is 
weak. We will typically consider jets that are denser than 
the lower surrounding medium (but also include models 
where the jet is already underdense in this region) . If the ex- 
ternal medium in the inner region is denser than the jet, one 
expects strong disruption of the jet and hence diffuse and 
destabilised jets in this inner (parsec scale) region where the 
jet would decelerate and dissipate its energy. This is in con- 
trast with observations of narrow and relativistic jets in the 
inner region of FR I galaxies and in the HYMORS, where 
the energy deposited at the la rge-scale is comparable to th e 
energy of the central engine (Rawlings & Saunders"1991|. 
In this work, we are particularly interested in radio loud 
FR I with a powerful jet at high Lorentz factor in the lower 
scale. This is the case in the HYMORS, where the jets have 
the same properties at small scale than the high energy FR 
II jets. The differences appear only in one side of the jet at 
the larger scale, where the jet morphology changes to show 
structure characteristic of an FR I. Therefore, we assume 
that further downstream, the jet encounters a high density 
medium, such that the jet undergoes a strong deceleration 
and a strong compression. This latter part is typically the 
only part simulated in previous numerical jet studies, where 
a prescribed jet configuration penetrates a usually uniform, 
high density, external medium. Moreover, we study the ef- 
fect of the initial opening angle of the jet in its interaction 
with (layered or stratified) external medium. 

For the density jump in the case of 3G 321, we can think 
of the lower region as the rarefied medium in the inter- 
cluster medium, while the upper region represents denser 
medium of the companion galaxy. In other HYMORS, the 
density jump is then thought to occur on one side of the 
AGN where there are denser molecular clouds in the inter- 
cluster environment or interstellar medium. One of the 
main points we hereby address for the first time is the 
change in jet head properties during the propagation phase 
in the inner region, and how this in turn affects the jet sta- 
bility in the upper, denser medium. Indeed, the jet inter- 
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acts with the inner medium mainly through the bow shock, 
whereas the jet beam is only weakly disturbed laterally. 
The shocked swept-up matter during this phase, and the 
shocked part of the beam, all constitute a structured bow 
shock ahead of the beam, and both shocked regions form 
a new hot layer with lower Lorentz factor, characterised 
by lower Mach number. The interaction of this preformed, 
structured jet head with denser external medium in the 
outer region should produce a strong cocoon and backflow. 
This latter will also disturb the non-shocked jet beam as 
it penetrates the denser medium. This will increase the en- 
trainment of ambient material through velocity shear in- 



stabihties and decelerate the jet ( ,De Young||1993 Perucho 
fc Marti|[2007l ). 



2.3. Initial conditions 

Deducing the precise internal properties for an FR I jet and 
its environment from observations is a difficult task, and 
density contrasts in particular have partly been obtained on 
the basis of numerical studies. However, from generic prop- 
erties of jet propagation in various FR I galaxies, one can 
estimate the kinetic luminosity and the jet beam Lorentz 
factor. Then, the choice of jet and environment parame- 
ters is determined on input by the kinetic luminosity of 
the jet and the estimated jet propagation speed in the two 
media. The kinetic lum inosity of a typical powerful jet is 
-^iet.Kin/- 10;^^ergs/s jR awlings fc Saunders" 1991' balyl 
T995l|Carilli fc Barth el|[ra96 Wan et al. 2000: Drake 200^ 
Tavecchio et al.||2 004; Kino & Takahara 2004) and we used 
this observationally supported value for the simulations in- 
dicated as cases A, B, G, H, I, J (see Table [T]). We in- 
vestigate also the propagation of low energy jets and com- 
pare them with these more powerful jets. For more radio 
quiet galaxies, a lo wer energy jet with Ljot.Kin ~ 10^'^ ergs/s 
( [Allen et al.||2006 ) is deduced and this value is used for the 
simulations C, D, E, F. The integrated energ y flux over 
the beam cross s e ction is computed f rom (e.g. | Bicknell fc' 
Begelmanl (|1996|); |Mart7et al., ( ,1997| ; [Rosen et al., ( ,1999, ); 
Scheck et al.| ( [20021 ] 



^jot.Kin 



(7b /ib - 1) /Cb7b7ri?b^b , 



(1) 



where the subscript "b" indicates the jet beam, pb is its 
density, 7b is the Lorentz factor, Vh is the speed, pb^b = 
Pb+ f~L ^^'^ enthalpy, and _Rb is the jet radius. For the 
latter, when we assume a jet with opening a ngle of 9h = 3° 



at 1 p c, the observed value in Centaurus A (Horiuchi et al. 
2006[ |, we get a jet radius Rh O.OSpc. We will tix the 
radius i?b = O.OSpc on our jet input boundary, which in 
turn is assumed to be located at a distance from the source 
of 0.5pc. 

Further, an estimated propagation speed of the head 
of the jet can be obtained from using the expression for 
pressure-matched jet propagation in ID ( Marti et al.|1997 



Rosen et al.||1999" I . For a cold external medium, this yields 



Vh , 



(2) 



where rj^i = 7b p^h" ratio between the inertia in 

the jet and in the external medium. In the relativistic case, 
the inertia of the flow increases as 7^ with the speed (which 
makes relativistic jets more stable than jets in young stellar 



objects). Our t = conditions then use the number density 
of the cold external medium in the pc-scale region set to 
JT-Low = lcm~^ as a scaling value, together with i?b as a unit 
of length, in combination with c ~ I. Then, the pressure in 
this cold external medium is set to p — 10~^TOp riLowC^, and 
this value for the pressure is in fact taken equal in the jet, 
and also in the outer region which is in static equilibrium 
with the inner region. We take the beam Lorentz factor 
fixed at 7b = 20 for the simulations A, B, E, F, J and at 
7b = 10 for the simulations C, D, G, H, L These assumed 
high inlet Lorentz factors are at the observed v alues appro- 
priate for FR II and BL Lac ob jects at pc-scale (Kellermann| 

This is consistent with tEe 
R I radio galaxies observed 



] [2004l ICoheri et al.|[2007| 
thesis that BL Lac axe F 



et al. 1 12004 

hypot 

with a small angle to the line of sight ( Urry fc Padovani 



1995) 



The density in the jet, and in principle also the pressure, 
can be deduced using Eq. [1] and Eq. [2] by imposing the 
jet head propagation Lorentz factor 7head,Low in the low 
density medium (in practice, this is limited by the condition 
to have a positive jet pressure and density). If we choose 
a value 7hoad,Low = 5 (i-head.Low = 0.979796), we find from 
Eq. [2[ a ratio between the jet beam inertia and the lower 
medium inertia. This is in the range 77r,low — [2672.3, 3000] 
for our models with high energy (A, B, G, H, I, J), while we 
have ?7r,low = [1.8,2.7] for our models with lower energy 
(C, D, E, F). We anticipate from these inertia contrasts that 
the high energy jet in the lower region is fairly stable and 
will conserve its narrow structure, while the lower energy 
jets will already be disturbed by the external medium in the 
lower part. When we use Eq.[T]to obtain our computational 
value for the density ratio between the jet beam and the 
external medium in the lower scale, we find rib/'T-Low — 
6.681 for the simulations A, B, J; rib/nLow — 29.99 for 
the simulations G, H, I; nb/riLow — 2.738 x 10~^ for the 
simulations C, D; and finally nb/^Low — 6.82 x 10~^ for 
the simulations E, F. 

For the density of the external medium in the upper, 
downstream region, we can argue similarly by setting a 
head propagation Lorentz factor at kpc-scale. In this pa- 
per, we will investigate many cases for the upper medium 
conditions. In the two first cases (A, B), we choose a uni- 
form density medium. In all other models, beyond the jump 
in the density at the interface between the lower and higher 
scale region, the density decreases with distance fr om the 



source with a simp le power-law ^jump/V-R^ + Z'^ (Kaiser 



. Furthermore, in model A, we assume 



that the jet undergoes weak deceleration in this upper re- 
gion, where the Lorentz factor of the jet head drops to 
7hoad,Up = 1-5 (mildly relativistic). Then, the ratio be- 
tween the jet beam inertia and upper medium inertia is 
'7R.Up = 8.65, which increases the influence of the external 
medium on the jet. Again, similar reasoning makes it plau- 
sible to use a density ratio nup/'T-jot ~ 4.687 for the upper 
medium then. In model B, we consider a very dense upper 
medium, expected to induce a strong deceleration of the jet 
to 7hcad,Up = 1-02 (i^hcad.Up = 0.197), which corresponds to 
a sub-relativistic jet. The inertia ratio between the jet beam 
and the upper medium is in this case ??r,up 0.0604, and 
this very low value will lead to increase the jet-external 
medium interaction, and the growth of body and surface 
mode instabilities in the jet. We find that to end up with 
this low Lorentz factor in the upper high density medium, 
we need to assume a fairly extreme, high value of density 
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Fig. 1. A zoom on the jet head at t = 160 light crossing times of the jet radius Rh, for cases A, B and J, just prior to 
penetrating the denser upper medium. Shown are: (Top) the logarithm of density, (middle) the Lorentz factor, (bottom) 
the effective polytropic index. In this and following figures, distances indicated on the axes are shown in parsec. 



'jet 



671.22. This is admittedly very high, however. 



many numerical simulations show that to decelerate a rel- 
ativistic jet to sub-relativistic speeds, we need overdense 
external medium where ^up/fib > 100 and even higher 



boundary 3°, since the jet is supposed to coUimate slowly 
during its propagation. 



( Krause 



20051. This is also seen by the exploited valued 2.4. Employed resolution 



for the density contrast in most numerical simulations to 
date, as collected in Table [T] Furthermore, from Eq. [TJone 
can argue that for fixed luminosity, the jet density decreases 
with jet radius as (0.05pc/i?b)^- This suggests that for big- 
ger opening angle of the jet (hence bigger jet radius), the jet 
density and thus also the upper external medium density 
will be lower, while the behaviour of the jet in this region 
should be the same. For all other models the density ratio 
at the interface reaches values as indicated in Table [l] and 
then decreases outwards. To model also the effect of conical 
versus cylindrical jet propagation, we assume a 1° opening 
angle at the inlet for models D, F, H, I and J. This opening 
angle is taken smaller than the open angle of the jet at the 



In all our simulations of jet propagation, we set the lower 
boundary at Zi^ = 0.5 (^g-^l^^ pc and initialize the jet ma- 
terial within the domain for radii i? < i?b and extended to 
Z — 1.0 ^ ?5pc ) P*^- ^^'^ boundary imposes a stationary 
mass flux at the lower boundary for radii R < Rh- The 
lower boundary for R > i?b is open. The top boundary is 
set at Zcxt and the jump in density is set at a distance 
Zjump — -^0x1/2, midway in our (very large) computational 
domain. Note again that, as evident from Table [T] we here 
cover jet propagation over distances that go as far as 400 
jet beam radii and with enough resolution to analyse with 
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high accuracy the interaction of the jet with the external 
medium, and this is only feasible thanks to our AMR ca- 
pabilities. 

Simulation A is done on a domain with size [R, Z] G 
[0, 10] X [10,400] (in units of i?b), with a resolution on the 
base level of 40 x 1200. Our grid-adaptive runs allow for 5 
levels, achieving an effective resolution of 640 x 19200. The 
second simulation is done on a domain with size [R, Z] e 
[0,40] X [10,400] and base level resolution of 160 x 1200, 
however, we allow for 6 levels achieving an effective reso- 
lution of 57 X 38400. This anticipates that since the ratio 
between the jet beam density and upper medium density is 
very high in the second case, the shock in this model will 
be very strong, and is likely to produce a turbulent cocoon 
and jet in this region. All other cases are done on a domain 
with size [R, Z] g [0, 40] x [10, 400] and with an effective res- 
olution 3072 X 4800 (4 levels). All the simulations exploit 
the hybrid version of HLLC, as explained in the appendix. 



3. Discussion of results 

From a basic point of view, the jet-external medium in- 
teraction is structured along the jet propagation axis in a 
way which is similar to the ID tests described in our ap- 
pendix. At the head of the jet, there are the forward shock, 
the contact discontinuity (called the work surface) , and the 
reverse shock (called the Mach disk). The forward shock 
compresses and heats the external medium, which spreads 
laterally, leading to the formation of a bow shock. The re- 
verse shock decelerates the beam matter by converting its 
kinetic energy to thermal energy. However, the shapes of the 
Mach disk and forward shock in a true 2D jet are oblique, 
and other 2D effects appear which we describe in the fol- 
lowing. 

3.1. Propagation through the uniform lower region 
3.1.1. Models A, B and J 

In models A, B, and J, the properties of the jet and of 
the lower external medium are the same, and the following 
discussion is applicable to these 3 cases. The jet density, 
prior to penetrating the dense upper region, is seen in a 
close-up view in Fig. [Tj This figure is at time t = 160 (in 
units of light crossing time for the jet beam radius). In the 
inner low density medium, the forward shock is relativistic. 
The jet also interacts laterally with the external medium, 
as there is a boundary shear layer. Due to the favorable 
density contrast, this layer does not disturb the jet propa- 
gation in this lower region. This thin shear layer is created 
as the Mach disk compresses the external shell of the jet. 
Only near the jet head, in between working surface and re- 
verse shock (Mach disk), and a bit beyond the reverse shock 
location, is the shear layer Kelvin-Helmholtz unstable, see 
our (zoomed-in) Fig. [l] We find hardly any backflow from 
the work surface in this lower region. 

Focusing on the axial structure of the jet head. Fig. [2] 
compares the obtained on-axis structure with a first equiva- 
lent ID Riemann problem. The sound speed in the shocked 
external medium (between forward shock and contact or 
work surface) increases to reach the maximum value al- 
lowed by the equation of state Cg ~ 0.636, hence the effec- 
tive polytropic index in this flow drops to Fes ~ 1.34 (a 
near ultra-relativistic state). Then the compression is very 
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Fig. 2. Top: A cut along the Z axis through the jet head at 
t = 160, prior to penetrating the denser medium (cases A, 
B and J). Shown is: (Left) the logarithm of density, (Right) 
the effective polytropic index. Bottom: a ID equivalent rel- 
ativistic shock tube problem. 



high and hence the distance between the contact disconti- 
nuity and forward shock remains small: we obtain a spacing 
of about 6X = 0.044 pc at t = 160. This is even three times 
smaller than that obtained in an equivalent ID shock prob- 
lem, as seen in Fig. [2j The difference is induced by the 
lateral spreading of the shocked shell of swept-up matter, 
which acts to decrease the amount of accreted mass. Indeed, 
shocked material in front of the beam will start to spread 
laterally, due to the strong (radial) gradient in the pres- 
sure and in the inertia between the shocked material and 
the cold, static external medium. This initial spreading can 
even be quantified from a second ID Riemann problem in 
the direction across the beam (which is not shown here), 
and this confirms the 2D result which shows an initial su- 
personic lateral speed wiatorai ^ 0.7973 for the shocked ma- 
terial. This outwards spreading leads to the formation of a 
bow shock, whose lateral spreading eventually (a 2D effect) 
decelerates to uiatorai ~ 0.36. With this fast lateral spread- 
ing of the shocked material, no complex cocoon forms in this 
region of low density for these three cases A, B and J. The 
second equivalent Riemann problem in the radial direction, 
also shows how a rarefaction wave propagates laterally in- 
ward into the swept-up matter by the jet. This rarefaction 
wave moves slowly towards the axis at approximate speed 
Vi- — 0.058, and matter accumulates mainly near the jet 
axis. This produces relatively more deceleration of the jet 
near the axis, and hence to a radially structured jet head. 

Returning to the axial structure as shown in Fig. 2] in 
this lower density region at t = 160, the reverse shock (or 
Mach disk) is Newtonian with ^ 1.013. This in fact, is 
totally different from the equivalent ID problem, as clearly 
seen in effective polytropic index in Fig. |2] The Mach disk 
compresses the beam matter, increasing the density and de- 
celerating the flow to 7 ~ 18 (in the equivalent ID case the 
shocked jet beam should have 7 = 5). Furthermore, the 2D 
result shows a relatively larger region of the jet head which 
is formed by shocked beam matter (compare the ID cuts 
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Fig. 3. A zoom on the jet head after penetrating the denser medium, for case A. We zoom with R, Z G [0,0.5] x [9, 13] 
parsec, at t = 240. Top: logarithm of density, middle: Lorentz factor, bottom: effective polytropic index. 



in Fig. |2]). The sound speed of the shocked beam material 
reaches ~ 0.2 with an effective polytropic index ~ 1.64. 
Also this shocked beam material tends to expand laterally, 
with a speed uiatorai '^-^ 0.1. Similar to the lateral dynamics 
of the swept-up, shocked matter, again a rarefaction wave 
propagates towards the axis, which progresses very slowly, 
and eventually the shocked jet head ends up denser and 
colder. Since the shocked external medium (in the forward 
bow shock) has already cleared the immediate surroundings 
of the beam, the shocked beam matter will hardly form a 
backfiow. 

The remarkable feature of jet beam propagation in a 
low density region is thus its stability. There is a very weak 
influence of the external medium on the jet, which allows 
the jet to transport energy over a long distance, and we 
followed this process over about 200 light crossing times of 
the jet beam radius. Only a small fraction of the jet energy 
is transferred to the external medium in the inner region, 
through the bow shock. During the propagation in the inner 



region, there is no strong backflow and we hardly see the 
formation of a turbulent shear layer. In this lower region, 
the Lorentz factor of the jet beam, and even the shocked jet 
beam, remains relativistic with 7 ~ 20. In turn, the jet head 
propagates relativistically with a Lorentz factor 5. It must 
be clear from the above discussion that even in the lower 
density region, the inclusion of a realistic EOS is necessary 
as the forward shock is relativistic and the reverse shock 
(Mach disk) is Newtonian. 



3.1.2. Other models when traversing the lower medium 

In model C, which is a light jet in this lower region at 
beam Lorentz factor 10, the jet interacts with the external 
medium mainly through the bow shock. Due to its lighter 
density, a turbulent cocoon gets formed during propagation. 
The turbulent cocoon interacts also laterally with the jet 
beam by compressing it, inducing internal shocks in the jet 
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beam. The distance between shocks is about Ashock — 20i?b- 
At each shock the jet undergoes a small deceleration to 
Lorentz factor 7 ~ 9 and is accelerated again to Lorentz 
factor 7 10 due to the expansion of the jet after each 
shock. Thus, in this lower region the jet remains relati vis- 
tic, its energy is mainly kinetic, and only a very weak frac- 
tion of the energy is transferred to the external medium 
by entrainment. The interaction of this light jet with the 
external medium also shows the formation of a low speed, 
higher pressure sheet of thickness ~ -Rb/2 surrounding the 
inner jet. 

Model D is similar to C, but has the jet emerging with 
an initial open angle of 6* = 1°. In the first part the D jet 
undergoes a free expansion, and its density and pressure 
decrease until the pressure in the surrounding cocoon be- 
comes high enough to compensate the jet thermal and ram 
pressure in the lateral direction. One can witness that in 
this region, an oblique shock forms, which moves with a 
slow speed w ~ 0.3 that decreases during the jet propaga- 
tion. In fact, as the jet propagates forward, the pressure 
of the cocoon behind the jet head decreases with cocoon 
expansion, and then the oblique shock moves forward. At 
t = 820, the oblique shock is at normalise distance Z 5 



Fig.(ll) This oblique shock in fact decelerates the jet to 
7b ~ 8 and consecutively coUimates it cylindrically, to end 
up with a smaller radius than the initial one imposed at in- 
let. Beyond the oblique shock, the jet beam is unstable and 
subject to the development of internal shocks. The amount 
of the matter entrained by the jet is larger than in the C 
case (cylindrical inlet). In this part of the jet, the inter- 
nal energy increases and the state of the matter becomes 
mildly relativistic with an effective polytropic index falling 
to Toff ~ 1.55. The jet in this model D is then more stable 
than the jet in the model C, since the initial conical propa- 
gation of the jet decreases the lateral interaction of the jet 
with the external medium and then the jet propagates in a 
laminar way until it reaches the coUimation shock. 

The model E differs from C only in its twice higher 
beam Lorentz factor. It has an inertia ratio between jet- 
external medium of similar order than in model C. Despite 
the faster jet beam, its Lorentz factor along the beam in 
the first region ends up oscillating between 10 — 25 due to 
the successive internal shocks and jet expansions, such that 
the jets in both model have the same general behavior. In 
the model F, the E (cylindrical) jet model now emerges 
conically with an initial open angle oi 9 = 1°. Its Lorentz 
factor 7b — 20 and kinetic luminosity ijct.Kin — 10'*'^ergs/s, 
makes that the difference with the conical model D is that 
the jet is faster. In this model F, the forward shock that 
forms at the jet head is stronger producing higher pressure 
and a more extended cocoon. This cocoon limits the region 
of the free conical expansion of the jet by forming oblique 
shocks that recoUimate the jet. This oblique shock in model 
F propagates with a speed v ~ 0.1, which is lower than the 
speed of the oblique shock in the model D. As in model D, 
the speed of the oblique shock decreases during the propa- 
gation. Beyond this oblique shock, the Lorentz factor of the 
jet decreases to about 10, and the jet starts to be unstable 
with the formation of internal shocks that compress and 
induce consecutive acceleration and deceleration of the jet. 
In this region the state of the matter in the jet is relativistic 
with an effective polytropic index Fefi ~ 1.45. At t — 820, 



In model G, the interaction between the jet and the 
external medium is roughly the same than in model A, be- 
cause only the Lorentz factor of the jet beam changes from 
7 = 20 in model A to 7 = 10. In fact, the main differences 
between the two models is in the details of the structure 
of the jet head, since the rate of compression at the front 
shock and the Mach disk (reverse shock) is different. In any 
case, the jet interacts weakly laterally with the external 
medium and it remains stable, consistent with its higher 
density with respect to the surroundings. Finally, in the 
conical models H and I (otherwise similar to G and also of 
high energy), the jet interacts with external medium only 
through the shock in front and there is no lateral interac- 
tion between the jet and the external medium, as only the 
head of the jet is a bit coUimated. In fact, the conical prop- 
agation of the jet decreases the influence of the low density 
external medium on the jet. 

3.1.3. Summary for lower region 

In the lower region, the jet head propagates with a speed 
Vhead.Low ~ 0.6 for the models C, D, E, F, and with a 



speed w/i, 



0.99 in the models A, B, G, H, I, J, 



the oblique shock is at normalise distance Z ~ 7.5 Fig.( 11 ) 



as collected also in Fig. [T2] In fact, these two groups are 
clearly distinct from the (input) kinetic luminosity of the 
jet beam. The low energy first group finds the lighter den- 
sity jets interacting strongly with the external medium and 
slowing down, since in this group the power of the jet beam 
is ^jct,Kin ~ lO^'^ergs/s and then the inertia ratio between 
the jet and the external medium is low, and at most of order 
r^R ~ 2.73. For the second high energy group, the external 
medium influences weakly the jet, since the power of the 
jet beam is ijct,Kin ~ 10^^ ergs/s and then the inertia ratio 
between the jet and the external medium is high, up to or- 
der ~ 3000. Thus the variation in the Lorentz factor 10 
to 20 of the jet beam and the small variation in the open- 
ing angle of the jet (cylindrical case versus conical with a 
small opening angle 1°) do not significantly influence the 
speed of propagation of the jet head in the lower scale re- 
gion, which is here simulated up to 200 jet radii. While we 
pointed out various aspects that are clearly captured only 
using a relatistic EOS model, our findings for jet propaga- 
tion through the uniform, lower region are fully consistent 
with earlier works. 



3.2. Upper region propagation 
3.2.1. Model A 

We now turn to the second stage in the dynamics, after 
the jet has passed the density jump. In the model A, the 
density ratio between the jet and upper external medium 
suddenly changes to pup/Pb ^ 4.687, making it a light jet. 
In this denser upper region, the initial interaction between 
jet and dense upper medium is shown in a zoomed Fig. |3] 
at time t = 240. First, as the jet penetrates the denser 
medium, an oblique shock propagates laterally in the upper 
medium. The relatively higher density of the upper medium 
gives rise to reflection of this shocked matter and hence 
produces a somewhat more turbulent and hot cocoon. At 
the jet surface, a thin, rarefied, and very hot region develops 
(seen best in the low effective polytropic index in Fig. pi). 
In fact, the head of the jet now sweeps up more matter, 
enhancing the temperature. The high pressure in this region 
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in turn gives rise to another oblique shock, which is weaker 
but propagates upstream and toward the axis. This confines 
the jet. This shock is located at time t = 240 at about 
Z ~ 9.8. The jet radius drops due to this compression to 
about R = 0.03. 

Again we can learn from the analogous ID Riemann 
problem and compare it with the axial structure of the 2D 
result. This is done in Fig. [4] at a time corresponding to 
Fig.|3] As a ID shock-structured jet head penetrates in high 
density medium in an upper region, there are in fact four 
layers of shocked material that can be distinguished along 
the beam axis. Once the previously formed forward (bow) 
shock meets up with the density jump, a new forward shock 
develops which seperates swept-up high density external 
medium from static external medium. A second contact or 
work surface from then on seperates shocked high density 
matter with previously swept-up, shocked, lower density 
matter. From this same location, a new reverse shock forms, 
which traverses the previously formed structured jet head 
(i.e. consisting of shocked lower matter, shocked beam mat- 
ter, and ending with the old Mach disk or reverse shock). 
These add up to 5 discontinuities in total, which are clearly 
seen in the analogous ID problem shown in Fig. [4] However, 
the 2D jet propagation is rather different from this ID 
model. The difference is quantified in Fig. [4] as we draw 
the cut along the Z axis (middle row), and also at a fixed 
radial distance away from the axis at i? = 0.02 (bottom). 
The latter radius still remains in the jet 'spine', since the 
laterally bounding shear region which we also can detect in 
our jet beam variation extends from 0.025 to 0.03 parsec. 

We can clearly see that, when comparing these three 
cases, in 2D we find faster beam flows immediately behind 
the front shock, which are also relativistic with a Lorentz 
factor 7^4. This is accompanied by strong compression, 
as the effective polytropic index drops to Fes ~ 1.37. This 
is as in the ID case. To be precise, at time t = 240, the front 
shock reaches in the 2D case a distance Z ~ 12.75pc, where 
the ID reaches Z ^ 12.8. Behind the front shock in the 2D 
case, various internal shocks develop, partly induced by the 
initial structure of the jet beam with its interaction with 
the denser medium. At each new oblique shock, the flow 
undergoes a weak acceleration behind it. In the region of 
the shocked beam, the Lorentz factor along the axis is oscil- 
lating between 7inin ~ 16 and 7max ~ 18. Now, the trailing 
reverse shock (Mach disk) and other shocks in front of it 
are Newtonian along the axis, and the effective polytropic 
index oscillates between Feff^min 1-65 and Feff^max ^ 1.61. 
In contrast, along the spine at fixed radius R = 0.02, the 
reverse shock (Mach disk) and all shocks in front of it are 
stronger. In fact, the Mach disk is near-Newtonian (effec- 
tive polytropic index Foff ^ 1.58), whereas the other shocks 
are relativistic with effective polytropic index oscillating be- 
tween 7i„in ~ 1.34 and 7max ~ 1-46. 

There are thus clear 2D effects, such as the lateral com- 
pression of the jet by the hot cocoon, and the presence of 
the oblique Mach disk. The latter at the same time shocks 
the jet beam, while compressing it laterally. A kind of re- 
flection happens (see also Fig. l3l, and in front of the re- 
flection point, the jet radius increases again, and gives rise 
to a weak acceleration. We find almost no backflow along 
the jet in this upper region. In fact, parallel to the jet, the 
hot shell of shocked external medium (as discussed above), 
moves in the same direction with a speed ~ 0.4. The 
difference of speed between this shell and the surrounding 




lo IS 14, 18 la ao 

z 



so — 
15 r 
lO r 

5 r 




lO IS 14, 18 la SO 

z 



Fig. 6. A cut through the jet for case A along the Z axis, 
a,t t = 380. (Top) the logarithm of density, (Bottom) the 
Lorentz factor. 



slower and denser cocoon, makes this region subject to in- 
stability, as seen clearly in Fig. [3] Meanwhile, the jet itself 
remains relatively stable. It seems that the formation of 
dense shear in the outer region of the jet increases the jet 
stability. However, in the jet, weak internal shocks are in- 
duced by the compression. The end result is the formation 
of a knot with long wavelength, and this result is shown 
in Fig. [5] where we plot the logarithm of the density and 
Lorentz factor at a later time t = 380. We also show the 
variation of density and Lorentz factor along a larger sec- 
tion of the jet axis at this endtime, in Fig. [6j It can be seen 
that several internal shocks have developed, reminiscent of 
knots. 

From this first model, it appears that a mildly overdense 
external medium can not induce the strong deceleration 
observed in FR I jets, as the jet head for this case A remains 
relativistic in this region. However, this model could be 
relevant for FR II jets, as it reproduces the weak interaction 
with the external medium in the inner region, and stays 
stably structured in the upper region. 

3.2.2. Upper region: model B 

In the model B, the density ratio between the jet and the 
denser external medium is increased to Pup/Pb = 671.220. 
The interaction of the jet with this very dense external 
medium is very strong and entails the formation of a tur- 
bulent cocoon. In Fig. [7] we show the density in the entire 
simulation region, at our endtime t = 900. In this model, 
the inertia ratio between jet and external medium was small 
T] ^ 0.0604, making the jet beam subject to instabilities. A 
zoom of the jet head at this same endtime is also shown in 
the last panel from Fig. |8] 

The result of the interaction is the formation of a bow 
shock propagating in front of the jet beam and to the 
side. The shocked external medium initially spreads lat- 
erally with a speed ?;iatcrai ~ 0.18, lower than the sound 
speed. The high density of the external medium leads to 
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Fig. 4. To analyze case A in the upper medium: we show at Left: logarithm of density, center: effective polytropic index, 
right: the Lorentz factor, at time t = 240, i.e. after penetrating the denser medium. The three rows correspond to: (Top) a 
ID equivalent relativistic shock problem, (middle) a cut along the Z axis (i.e. R—0) through the jet head, and (Bottom) 
a cut along the axis at a radius R = 0.02. From the top to the bottom, note the scales difference: Z £ [12.0, 13.0] for top 
and middle, while bottom panel has Z G [11.0, 13.0]. 



a weak compression rate, and in fact, the front (forward) 
shock is now Newtonian with an effective polytropic index 
Foff ^ 1.6. The compression rate is about 5 and the jet head 
propagates with a Lorentz factor less than 1.25 along the 
axis and more slowly away from the axis. Indeed, along the 
line R — 0.01 the shock propagates with a Lorentz factor 
1.09 (Fig. [9]). 

The slow spreading of shocked dense matter makes the 
2D effects on the propagation of the forward shock less sig- 
nificant than in the first model, and there is less difference 
with a similar ID case. Again, Fig. [9] compares analogous 
ID results with axial cuts along axis and at R = 0.01, for 
an earlier time t — 240 (for which 2D impressions are shown 
in the top two panels of Fig. Isl. The shocked jet beam in 
2D is structured axially as in the ID case with the appear- 
ance of a new work surface, and the Mach disk, while the 
compression rate at the trailing Mach disk approaches the 
ID case. However, we see that the effective polytropic in- 
dex there is higher, namely Fcff ~ 1.4, and the Lorentz 
factor is different (7 ^ 1.03 in the ID and 7 ~ 6 in the 2D 
jet). The resemblance of the structure of the shocked ex- 
ternal medium and the jet beam in ID versus 2D is better 
near the jet axis. We can clearly detect the new forward 
shock, and the second contact discontinuity, seperating the 
shocked high density matter with previously swept-up mat- 
ter from the lower density region. The differences between 
2D and ID are the additional oscillations on the density 
and pressure, which are due to the lateral effects. Sideways 



spreading matter is strongly reflected, as the high density of 
the external medium slows down the lateral growth of the 
cocoon. Moreover, the oblique Mach disk confines the jet. 
Beyond it, the jet radius increases and the flow is acceler- 
ated again to reach a Lorentz factor 7 ~ 17. Then, the hot 
and dense cocoon compresses the jet again, and a new inter- 
nal shock decelerates the jet. Along the line R — 0.01, the 
reverse shock is Newtonian and there is evidence of more 
internal shock development, indicating complex interaction 
with the cocoon. 

In this case B, we do find strong backflows. The shocked 
jet beam matter by the trailing Mach disk again is subject 
to sideways spreading. The denser hot cocoon reflects this 
spreading matter, and induces the formation of a back- 
flow. In fact, between the jet and the shocked external 
medium, a shear layer develops, consisting of a dense and 
relatively hot flow propagating upstream with a speed of 
order Wbackflow ~ —0.2. Near the top of the jet head, the 
backflow is essentialy consisting of shocked jet beam mat- 
ter, flowing in a cylindrical shell of thickness AR ~ 0.01 
parsec. The mass flux in this backflow increases in the 
upstream direction. Moreover, this backflow entrains with 
it shocked external cocoon matter, and then carries more 
mass, which slows down its speed to z'backflow ^ — 0.1 when 
it reacheas Z = 10.0 (the original position of the density 
jump in the external medium). The backflow interacts also 
strongly with the non-shocked jet beam, and it compresses 
the jet beam, and induces shocks. The final result of this 
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Fig. 5. A zoom on the jet head with R, Z £ [0,0.2] x [17.5,20], at t = 380, for case A. (Top) Logarithm of density, 
(middle) Lorentz factor, (bottom) effective polytropic index. 



fairly complex interaction is confinement and overall de- 
celeration of the jet in the upper region. The jet radius 
decreases to i? ~ 0.02 in Z - 5.0, at t - 900.0. Beyond this 
point the jet radii starts to increase. And at the interface 
between the two medium a strong oblique shock form. The 
jet is strongly decelerate at this shock to 7 ^ 2.0 and the 
it pressure increases, such the relativistic mach number fall 
to Af '--^ 5. The high pressure jet expand lateraly in this 
region and to radius R 6 x Rjct.h- 

The interaction of the backflow with the shocked ex- 
ternal medium in the cocoon induces the development of 
Kelvin-Helmholtz instabilities. This is the result of the ve- 
locity gradient between the backflow and the cocoon and 
also because the backflow is slightly underdcnse with re- 
spect to the cocoon. It is clear from Fig.|8]that the entire re- 
gion is highly structured due to instability development. In 
this upper region, the density ratio was pup/Pb = 671.220, 
and we indeed flnd a transition to a beam suddenly trav- 
eling at a Lorentz factor of 7up ~ 1.25. This value is a bit 
higher than the value used for setting up the initial condi- 
tion (where we used 1.02). This difference can be explained, 
since the head of the jet beam is made up of hot swept-up 
matter and hot shocked beam matter, while the estimate 
assumed cold conditions. Moreover, in this upper region, 
the forward bow shock is now Newtonian. The various lay- 
ers at the head of the jet give rise to strong turbulence 
development. As result of this interaction and the entrain- 
ment of the externa matter by the jet is the deceleration of 
the jet in the upper region to 7 ~ 1.5 (u ~ 0.3). 

When estimating the overall energy budget, we find that 
in the jet interaction with the upper medium, about of 58% 



jet energy gets deposited in the upper lobe, while a fraction 
of order 10% is reflected in the lower region. 

3.2.3. Other cases: effects of opening angle and density 
decrease 

Low energy jets. 

In the model C, the density ratio between the jet and 
the dense external medium is Pb/pup = 4.91 x 10"'^. There 
is strong interaction of this light jet with the external 
medium, and this implies the formation of a high pres- 
sure cocoon. This turbulent cocoon disturbs the jet by in- 
ducing a strong backflow. The difference in pressure be- 
tween turbulent cocoon and the colder lower region pro- 
duces also backflows propagating into the lower region. 
At the interface between the two regions, the jet beam is 
compressed. The jet becomes unstable and its pressure in- 
creases. However, the jet remains relativistic in the upper 
region, since it decelerates only to a Lorentz factor 7^8. 
The decrease of the density in the external medium as as- 
sumed for this jet (see Table [T]) and the (minor) decelera- 
tion of the jet makes the forward shock weaker and then 
the pressure of the cocoon at the head of the jet decreases. 
Therefore, the jet beam radius actually increases slowly 
during propagation in this upper stratified region. 

In the model D, the jet beam that finally emerges 
from the lower region is coUimated cylindrically (despite 
its original opening angle, see the discussion in the previ- 
ous paragraphs) and ends up with a relativistic Mach num- 
ber M ~ 20 and with a Lorentz factor of 7 8. This jet 
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Fig. 8. A zoom on the jet head after penetrating the high density medium in case B, for R, Z G [0, 0.2] x [9, 12] Top 
and middle panel show at < = 240: (Top) Logarithm of density, (middle) Lorentz factor. Bottom panel shows, later at 
t — 300, the effective polytropic index, zoomed with R, Z E [0, 0.3] x [10.5, 13.5]. 



interacts again strongly with the external, upper medium 
forming a high pressure cocoon that induces a backflow, 
also propagating upstream in the low density region. This 
backflow compresses the jet near the (perturbed) interface, 
increasing its density and pressure. The jet in the upper 
region is found to decelerate to a Lorentz factor 7 ~ 5 
through multiple internal shocks. As in model C, the com- 
pression of the jet at the interface increases its pressure 
and the assumed density decrease in the external medium 
causes the jet to undergo a slight conical expansion in this 
upper medium. 

In the model E (main difference with C is that the jet 
is faster), we find an even stronger interaction with the 
external medium. This configuration enhances the develop- 
ment of the shear instabilities which causes entrainment of 
ambient material. Then this jet E is found to decelerate 



smoothly to a Lorentz factor 7 ~ 5, while it was character- 
ized by 7 '--^ 20 on inlet. 

In model F, the jet beam emerges from the upstream, 
lower region with a fairly low mach number M ~ 15 
(lower than in model D) and with an internal energy 
of the order of the mass energy. Moreover, the jet in F 
propagates with a high Lorentz factor and higher thermal 
energy. When it interacts with the downstream, denser 
medium, it forms a flattened (compressed) and turbulent 
cocoon. Then a stronger backflow develops than in case 
D, that propagates upstream in the low density region. 
This backflow this time eventually induces a deceleration 
of the jet to a Lorentz factor 7 ~ 10 (along the axis). 
At the interface, the jet pressure ends up higher than 
the pressure of the cocoon, and also this jet expands 
in the upper region in a conical shape until it reaches 
pressure equilibrium with the cocoon in this region. The 
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Fig. 9. Shown are: (Left) logarithm of density, (center) effective polytropic index, (Right) Lorentz factor, for three cases. 
The three rows are (Top) a ID equivalent relativistic shock problem, (middle) a cut along the Z axis on the case B jet 
head at t — 240. (bottom) A cut along the axis at fixed R — 0.01. 



jet interacts strongly laterally with the cocoon in this 
region, and overall, the jet beam slows to a Lorentz factor 
7 5 (a long the axis) through multiple internal shocks 
induced by instabilities and entrainment of external matter. 

High energy jets. 

In model G, the difference with model A is that the 
jet is slower and the upper medium less dense. The inter- 
action of the G jet with the external medium is therefore 
weaker than in A, and the jet also remains relativistic. At 
the end of simulation, this jet deposited 2.5% of its energy 
in the upper medium and a fraction of the order of 1.9% 
is reflected in the lower region. In model H, which already 
had a conical expansion of the jet in the lower region, and 
when the jet reaches the interface between the two regions, 
the density of the external medium is 155.5 higher than 
the density of the jet. When the jet goes into the upper 
medium, a strong bow shock forms and a high pressure 
cocoon develops. We also find a backflow that reaches a 
speed around Wbackfiow ^ —0.5 at the location of the inter- 
face. A shock wave forms there and propagates upstream 
in the lower density region. This shock wave compresses 
the jet beam and coUimates it cylindrically near the inter- 
face between the two medium. At t = 820, the this shock 
wave reaches the normalise distance Z ~ 5.0 Fig.(10l. In 



the upper medium, the high pressure cocoon induces the 
formation of an oblique shock at the jet head. This com- 
presses and the jet radius changes from i? ~ 3.0 x i?jet,b 
to i? -Rjct,b/2.0. In all, the jet decelerates abruptly from 



Lorentz factor 7 ~ 10 to 7 ~ 2 at this location. At this 
shock, the pressure becomes higher than the pressure in 
the cocoon and the state of the matter becomes relativistic 
and its effective polytropic index falls to F ~ 1.34. Beyond 
this shock the jet beam spreads with an angle of about 
^ ~ 6° (mostly influenced by the decreasing density) and is 
accelerated to a Lorentz factor 7 ~ 8 before it reaches the 
bow shock. At t — 820, the oblique shock is at normalise 
distance Z ~ 11.0 Fig. (10 1 At the end of simulation, the 



H jet deposited 40% of its energy in the upper lobe and a 
fraction of order 16% is reflected in the lower region. 

In the model I, the difference with the model H is that 
the upper medium is less dense than the jet. In the ini- 
tial periode, the jet interacts with the external medium 
mainly through the front shock until the bow shock starts 
to develop prominently in the upper region and its pressure 
increases. The high pressure of the cocoon compresses the 
jet and induces an internal shock in the jet that collimates 
it. Beyond this shock, the jet becomes unstable and mul- 
tiple shocks develop, decelerating the jet to Lorentz factor 
7 ~ 8. At t = 480, the oblique shock is at normalise dis- 
tance Z ^ 14.0 Fig.( 10 1 At the end of simulation, this I jet 



deposited 18% of its energy in the upper lobe and a fraction 
of order 8% is reflected in the lower region. 

In model J, the difference with model H is that the 
jet is faster (Lorentz factor 20 on inlet) and the external 
upper medium has a density 2.5 lower than in model H. 
The initial phase of the interaction of the jet with the up- 
per external medium produces also a shock wave propa- 
gating upstream, coUimating the jet cylindrically. Similar 
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Fig. 7. The logarithm of density a,t t = 900 on the entire 
simulation domain, for the high density upper medium case 
B. 



Lorentz factor 7 ~ 10 through the shock. Beyond this, the 
jet spreads again and undergoes multiple internal shocks. 
Ultimately, this decelerates the jet beam to Lorentz factor 
7 ~ 5. At t = 480, the oblique shock is at normalise dis- 
tance Z ~ 12.5 Fig.(10l The faster jet J propagates for a 
longer distance than tEe jet in H. In fact, jet J traverses a 
longer distance in free ballistical propagation since the ram 
pressure in the jet beam is higher. Then the oblique shock 
forms farther away than in H and its overall compression 
rate is lower. At the end of simulation, the jet deposited 
about 18% of its energy in the upper lobe and a fraction of 
order 12% gets reflected. 



3.2.4. Summary of all cases 



Figures Tojll gives an impression of the jet morphologies, 
at the end of the simulations. In our models A, G with 
high energy and upper medium of about the jet density, we 
end up with a relativistic jet in both regions (inner/upper). 
In the upper, slightly overdense region (for A), the inter- 
action with the external medium is stronger than in the 
underdense inner region, and we find the formation of a 
relatively smooth cocoon. In the model B with high energy 
and a much higher density in the upper medium, the jet is 
relativistic in the lower region, where it is seen to propa- 
gate with a Lorentz factor of about 5, and turns to a sub- 
relativistic jet with propagation speed 0.3 in the very dense 
upper region. In all models with low energy C, E, the jet 
remains relativistic in the lower region despite fairly turbu- 
lent cocoons, and is strongly decelerated when it crosses the 
interface between the two media. In the cases undergoing 
strong deceleration, a shock wave propagates upstream (a 
reflection from impacting the higher density environment), 
and this compresses and decelerates the jet. In the model 
with an initial opening angle, an oblique shock forms be- 
fore the jet reaches the contact interface between the two 
media for the cases D, F with low energy. This shock sur- 
vives the interaction with, and gets enhanced by the shock 
wave propagating upstream that forms when the jet inter- 
acts with the upper medium. In conical models H, I, J with 
higher energy, such an oblique shock forms only when the 
jet starts to interact with the higher density upper medium. 
These pronounced differences in the propagation behavior 
between the two jets, is clearly seen when we quantify the 
propagation of the jet head, which is collected as function 
of time in Fig. 12 for all models. We can conclude that 
this is relevant to the FR I/FR II dichotomy: we find FR II 
behavior where the energy of the jet is high and the density 
of the external medium is not high enough to significantly 
impart a strong refiected shock, and FR I behavior when it 
can, even if the density beyond the interface decreases and 
also if the jet is not coUimated. 



to model H, the high pressure cocoon that develops in the 
upper region causes a backflow that propagates upstream. 
This disturbs the jet there and induces the formation of an 
oblique shock that compresses, coUimates, and decelerates 
the jet. In fact, at the shock the jet radius now falls from 
i? ~ 4.0 X i?jet,b to i? '--^ ^jct,b- The jet is decelerated to a 



Finally, we mention that in all these simulation, when 
the jet crosses the intersection between the two medium, we 
witness the appearance of Richtmyer-Meshkov instabilities 
at the shocked contact interface. We were able to resolve the 
local vortical development, using high resolution. In fact, 
the Richtmyer-Meshkov instabilities evolve out of vorticity 
deposition on the interface, which is due to the propagating 
jet bow shock passing it. 




Fig. 10. Contours of logarithm of density for high energy simulations: G at t=380, H at t=800, I at t— 480, J at t= 480. 
R and Z are normalised to (20 x i?b). 




Fig. 11. A contour of the logarithm of density for low energy simulations: C at t=820, D at t=820, E at t=820, F at 
t=820. R and Z are normalised to (20 x i?b)- 
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Fig. 12. The position of the jet head as function of time, 
for all cases. 



4. Conclusions 

In this paper, we extended, presented and appl ied the rela- 
tivistic hydrodynamics AMRVA C code (Mchani "e"t al.|2007 



van der Hoist fc Keppens|2007 1 with a generalized variable 



polytropic index equation of state for the purpose of mod- 
eling the relativistic shocks in GRBs and AGN jets. As 
follows in the appendix, we used various shock-capturing 
schemes on a variety of test cases for adiabatic and non- 
adiabatic cases for code validation. We demonstrated the 
code ability with stringent new test problems, developed 
according to the astrophysical context. We tested the code 
also for a case with a heated flow, using the Synge-like 
equation of state from [Meliani et al. ( 2004 1 . The Riemann 
problems were also solved exactly, and the AMR simula- 
tions handled cases with Lorentz factors of order 100 accu- 
rately. 

We explored a new scenario for sudden deceleration of 
relativistic jets in the FR I radio galaxies. This model for 
jet deceleration is based on a density jump in the exter- 
nal medium, with density suddenly increasing to the upper 
medium. We include models of conical jets with an ini- 
tial opening angle, and allow for decreasing denisty profiles 
within the upper medium. We investigated the propaga- 
tion and the dynamics of relativistic AGN jets over a long 
distance, resolving small scale instabilities which develop 
in the cocoon and which could be responsible for particle 
acceleration. This study was only possible with adaptive 
mesh refinement. We quantified and discussed the deceler- 
ation of the jet resulting from its interaction with a layered 
medium. We point out that jet deceleration can be signifi- 
cantly aided by an internal oblique shock that forms in jets 
with conical expansion. Under fairly extreme (but not as 
extreme as pursued in other simulations to date) density 
conditions in the dense upper region, it can reproduce the 
strong deceleration observed in FR I jets. 

In the early phase, the head of the cold fast (beam 
Lorentz factor 10 or 20) jet propagates in the low den- 
sity medium at a high Lorentz factor 7 = 5. In cases with 
high energy, almost no cocoon or backflows develop in this 
phase, and as the jet is heavier than the external region, 
it behaves more ballistically. When observed in this region, 
the relativistic beaming is high, so that one would observe 
a one-sided jet. In low energy jets, a cocoon always forms 
and slow backflows develop. They in turn induce multiple 
internal shocks making the jet decelerate and re-accelerate. 



In jets with low energy and with an initial conical flow, 
the flow remains ballistic until it reaches a self-consistently 
forming oblique shock. This compresses and decelerates the 
jet and two regions appear in the jet beam, the first with 
high relativistic beaming, and the second with low Lorentz 
factor, showing clear disturbances by the surrounding co- 
coon. 

In the outer region where the density increases suddenly, 
we find increased entrainment of ambient material through 
velocity shear instabilities at the jet boundary and working 
surfaces. Note that we here for the first time address how 
the prior interaction of the jet with the low density ambi- 
ent medium in the inner region (which makes the head of 
the jet consist of swept-up ambient medium and a shocked 
beam) affects jet propagation in denser media. The pre- 
structured jet head has lower Mach number, which gives rise 
to a strong interaction with the denser ambient medium. 

For models with a weak increase in density in the ex- 
ternal medium, the cocoon shows little turbulent structure 
and no backflow appears. The jet in this case remains stable 
and relativistic in the denser upper region. Also a fraction 
of its energy is deposited in the external medium, but only 
through the frontal shock. This model could be relevant 
for jets in FR II. With this scenario, one can explain how 
most of the energy of the jet is deposited in the outer re- 
gion, since relatively little jet energy gets transferred to the 
medium in the inner region. 

Those models with more extreme increases in density in 
the external medium, showed the formation of an overpres- 
sured and turbulent cocoon. A strong backflow develops in 
the outer region, which disturbs the jet structure. The back- 
flow compresses the jet and induces internal shocks which 
give rise to knot formation. The jet becomes subrelativis- 
tic within the simulated domain. The models with initial 
conical jets show the development of an oblique shock that 
decelerates the jet. In future work, we will follow these jets 
over increasingly larger distances in 3D scenarios. We also 
intend to provide synthetic observational radio maps from 
our computed jet models. 

Appendix A: Relativistic hydrodynamics and the 
equation of state 

The special relativistic hydrodynamic evolution of a per- 
fect fluid is governed by the conservation of the number of 
particles, and energy-momentum conservation. These two 
conservation laws can be written as 



(pu^) =0, (Tn, = 0. 



(A.l) 



where p, u — (7,71?), and T''^ — phu'^'-u'^ + pg^'''^ de- 
fine, respectively, the proper density, the four- velocity and 
the stress-energy tensor of the perfect fluid. Proper density 
is related to the number density n in the fluid rest frame 
p = mup, where rrip indicates the particle (proton) rest 
mass. The definitions involve the Lorentz factor 7, the fluid 
pressure p, and the relativistic specific enthalpy h. For the 
(inverse) metric g^'^, we take the Minkowski metric. Units 
are taken where the light speed equals unity. 

These equations can be written in conservative form 
involving the Cartesian coordinate axes and the time axis 
of a fixed 'lab' Lorentzian reference frame as 



dU 

Ik 



E 



&7 



0. 



(A.2) 



Z. Meliani et al.: FR I Jet deceleration 



19 



The conserved variables can be taken as 

U = [D = J p, S = hv,T ^ -f^p h~ p — 7p] , (A. 3) 

and the fluxes are then given by 

F= [p^Vjj'^phvv + plj'y'^phv — jpv\'^ , (A. 4) 

where I is the 3x3 identity matrix. To close this system 
of equations, we can use the Synge-l ike equation of stat e 
(EOS) for an ideal gas as also used by Meliani et al. (20041, 
which is a poly tropic equation with a corresponding classi- 
cal polytropic index T relating 
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(A.5) 



where e 



eth is the specific internal energy includ- 



ing rest mass, and eth is the specific thermal energy. This 
specific internal energy is given in function of the pressure 
and density by rewriting the above equation to 
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The relativistic specific enthalpy is then given by 
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The sound speed (in light speed units) is then found from 



1 dp 
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At each time step in the numerical integration, the prim- 
itive variables {p,v,p) involved in flux expressions should 
be derived from the conservative variables U resulting in a 
system of nonlinear equations. One can bring this system 
into a single equation for the pressure p, directly following 
from the deflnition of the conserved variable r from 



iT + p + D)il-v{p)^) ~ ph{p)^0, 



(A.9) 



which, once solved for p yields v = S/{t + p + D) . One 
inserts h from Eq. (A.7 1, in which one uses Eq. ( |A.6 l, with 
in addition D = Tpwhil e 1/7 ^ = 1 - S'^ /{t + p + D)^. 
This nonlinear equation (A.9 1 is solved using a Newton- 
Raphson algorithm. For a chosen fixed index T = 5/3, we 
then achieve a locally varying, effective polytropic index 
taking on values between its relativistic 4/3 and classical 
5/3 extremes, found from 



- off 



r - 1 
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(A.IO) 



This provides an excellent approximation to the true Synge 
gas expression, while avoiding costly Bessel function evalu- 
ations. If 1 < r < 5/3, we incorporate effects ranging from 
near-isothermal conditions to adiabatic flow, and if F > 5/3 
we incorporate the effect of energy losses. 



Appendix B: Numerical tests for realistic EOS 

The main aim of this appendix is to point out that we have 
determined the exact solutions to newly selected Riemann 
problems with a realistic equation state. This exact solution 
is here used, and can be used by other authors, to bench- 
mark codes. We test the code for adiabatic EOS Riemann 
proble ms, as also investigated by Mi gnone fc McKinney| 
(2007), and for a more general case with equivalent classical 
polyt ropic index 3/2 to mimic the heating in a relativistic 
fluid (Meliani et al. 2004 1. In each test, we analyse the effect 
of the equation of state with variable effective polytropic in- 
dex on the behaviour of the shock and rarefaction waves. 
We include tests with high Forentz factor 100, with the 
aim to investigate shocks in Gamma Ray Bursts, where the 
forward shock is always relativistic and the reverse shock 
Newtonian. The matter shocked by the forward shock has 
then an effective polytropic index reaching 4/3. However, 
the matter shocked by the reverse shock has effective poly- 
tropic index 5/3. This difference state of the matter between 
the two shocked media induces a change in the propagation 
speed of the shocks. In fact, using this EOS will improve 
the simulations for GRBs and induce some variation from 



(A.6) the results iniMeliani et al.|(|2007l. 



B.l. Exact solution to Riemann problems with realistic EOS 

The solution of the one-dimensional Riemann problem in 
hydrodynamics consists of determining the temporal evo- 
lution of a fluid which, at some initial time, has two ad- 
jacent uniform states characterized by different values of 
uniform velocity, pressure and density. These initial condi- 
tions completely determine the way in which the disconti- 
nuity will decay after removal of the barrier separating the 
initial "left" and "right" states. 

In general, the Riemann problem requires the solution 
of a nonlinear algebraic system of equations written as a 
function of a set of unknown quantities. In relativistic hy- 
drodynamics (RHD) an exact solution has been obtained 
only rather recently and was proposed by |Marti fc Miiller] 
(1994^ for flows that are purely along the direction normal 
to the initial discontinuity. This work has then been ex- 
tended to th e case in which tangential velocities are present 
(Pons et al. ( j2000) ) and improved in efficiency by exploit- 
ing the relativistic invariant relative velocity bet ween the 



two states to pre dict the wave pattern prod uced ( Rezzolla 
fc Zanotti] ([2001[) and [Rezzolla et al.|p003l)). 



The exact solution is found after expressing all of the 
quantities behind each wave as functions of the value of the 
unknown gas pressure p at the contact discontinuity. In this 
way, the problem is reduced to the search for the value of 
the pressure that satisfles the jump conditions at the con- 
tact discontinuity. Recently |Giacomazzo fc Rezzo lla (20061 
found for the first time (after a first attempt by [Romero et] 



al. (2005 ) limited to a very special case) a general procedure 
to compute the exact solution of the Riemann problem in 
relativistic magnetohydrodynamics (RMHD). This solver 
can also be used in RHD by setting to zero the magnetic 
fiel d and in thi s case it is similar to the method developed 
by Pons et al] (|2b00). Moreover, the exact solver can im- 
plement, both in RHD and in RMHD, various equations 
of state of the form p = p(p, eth) making it a very general 
tool that has been used for the testing of several special 
and general relativistic numerical codes. To the best of our 
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Fig. B.l. One-dimensional Riemann problem for colliding, identical cold flows, with a jump in tangential velocity v., 
The solid lines are the exact solution. The result is plotted at t = 2. 




Fig. B.2. One-dimensional Riemann problem for two states with different initial thermodynamic state. The fluid at left 

has a classical state and the fluid at right a relativistic state. The solid lines are the exact solution. The result is plotted 
at t = 1. 



knowledge, this is the first time in which exact solutions 

with generic initial states are computed with a different j 

than polytropic EOS. I 

1.60 r 



5 1.50E 

B.2. Numerical schemes , , 



1.30 ; 



We include in this section various stringent tests to vali- 
date the code. We perform one-dimensional test problems 
to show the code ability to resolve Riemann problems, and 
show that we reproduce the exact wave patterns in spe- 
cial relativistic hydro with varying EOS. We use different 
discretization s chemes, namely TVDLF (Toth fc Odstrcil 



1996) , HLLE ( Harten, Lax fc van Leer| 1983| [Einfeldt 
1988|, HLLC (Mignone & Bodo 2005), andaEo a hybrid 



mixture between TVDLF and HLLC. The latter ensures 
that in regions where spurious oscillations can be induced 
by HLLC, the code switches to TVDLF. This switch is ef- 
fective in cells where the fluxes as computed at left and 
right edge change direction. For the linear reconstruction 
from cell center to edge, we use a robust minmod scheme 
that reduces any spurious numerical oscillations. All these 
ID tests are compared to their exact solution, computed as 
outlined in Section IB. II 



-0.3 



-0.1 



0.0 

X 



0.1 



0.2 



Fig. B.3. Riemann problem for two states with different 
initial thermodynamic state, with a stationary contact dis- 
continuity at X = 0. 



B.3. ID Tests: Riemann problems 

We start with the collision of two cold flows with, in addi- 
tion, an opposite orientation of tangential velocity. The ini- 
tial states are pL = 100, pl = 1, Wx, l = 5/-\/26, Uy^L = 0.01 
(left) and pR = 100, jjr = 1, Wx,r = Vy,R = 

—0.01 (right). Both flows have a Lorentz factor 7 '--^ 5 and 
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Fig. B.4. Onc-dimcnsional Ricmann problem for the inter- 
nal shock in a heated wind, colliding with slower outflow, 
plotted at i = 0.4. 



the classical polytropic index 5/3. The simulation is done 
with 20 cells at the lowest grid level, and we allow for 7 
levels on the spatial range —2 < x <2. The test shows the 
characteristic pattern of two shock waves and a station- 
ary discontinuity in the tangential velocity. The two shocks 
convert the kinetic energy to thermal energy and the ef- 
fective polytropic index which was Toff = 1.657 through- 
out drops locally to 1.345 increasing the compression rate 
Pel Pi. ~ 20. 44, w here pc denotes the central proper den- 
sity. In Fig. |B.1| we plot the result dX t — 2 as obtained 
with the hybrid version of HLLC. This scheme captures the 
two shock waves accurately, as they are propagating with a 
shock speed Wsh ~ 0.33, and seperate the two regions with 
different states of the matter. Indeed, a classical fluid exists 
in front of the shock waves and an ultra-relativistic state is 
found between the shock waves. 

In the second test, we follow the evolution of an ini- 
tial discontinuity between two fluids with different thermo- 
dynamic states, both moving to the right. At right, the 
relativistic fast gas is characterized by = 5000, — 
0.01, ■i;x,R = 0.5, = 0.5, and this interacts with the cold 
slower gas at left pl = 0.01, pl = 1, t'x,L = 0.1, Vy^i^ = 0.5, 
thereby creating a flow to the left. We take T — 5/3, and 
obtain at the left an effective polytropic index Fcfrx = 5/3 
and at right, it is Feff^R = 4/3. The simulation is done with 
20 cells at the lowest grid level, and we allow for 9 levels on 
the spatial range — 2 < a; < 2. The evolution leads to the 
formation of a shock wave propagating to the left, which 
heats and compresses the fluid. Between this shock and the 
contact discontinuity, the effective polytropic index of the 
compressed fluid drops to 1.45. In Fig. |B.2| we plot the re- 
sult at t = 1 as obtained with the hybrid version of HLLC. 
The compression rate pl,s/pl ~ 6.69, where pL,s denotes 
the density of the compressed fluid. More to the right of the 
contact, a rarefaction wave propagates into the relativistic 
fast fluid. A numerical difficulty is produced by the strong 
drop in the mass flux, even though the tangential velocity 
increases at the contact discontinuity, such that the Lorentz 
factor between the contact discontinuity and the tail of the 
rarefaction wave is 4.6 3 hig her than in front of the contact 
discontinuity (see Fig. B.2|. In this test, it is vital to intro- 
duce a realistic equation of state with effective polytropic 
index becoming function of the temperature. As the shock 
propagates to the left, it produces a mildly-relativistic adia- 
batic compression of the flow, where the state of the fluid is 



described by 4/3 < F < 5/3. This state can not be analysed 
with a classical constant polytropic equation of state. 

The third test represents an isolated contact disconti- 
nuity with only a jump in the density p such that there 
is also a jump in the effective polytropic index. The initial 
states are pL = 10^, pl = 1-0, Wx,l = 0.0, = 0.4 (left) 
and PR = 0.125, pr = 1.0, v^ji = 0.0, Vyji = 0.4 (right). 
The simulation is done with 12 cells at the lowest grid level, 
and we a llow for 2 levels on spatial grid —0.2 < x < 0.2. In 
Fig. B.3 we plot the result (at t = 2, representative of all 
times) as obtained with the HLLC scheme. The test shows 
that HLLC resolves an isolated stationary discontinuity ex- 
actly (and thus the isolated jump in effective polytropic 
index) . 

In the fourth test, we follow the evolution of an initial 
discontinuity between two flows to the right, with different 
Lorentz factor, both having a polytropic index F = 3/2. 
This represents the critical v alue in the Parker model for 
the solar wind (Parker 19601, and due to our EOS, both 
flows undergo a heating. The initial states are pL = 10.0, 
PL = 100.0, Wx,L = 0.866025, Uy,L = Cs,l (left) and 
PR = 1.0, PR = 10-3, i;x,R = 0.30491, z;y,R = C,,r (right) 
where Cg denotes the sound speed. The Lorentz factor at 
left is already order 100 here. The aim of this test is to mim- 
ick the internal shock between a mildly relativistic outflow 
and subsequent hot and fast ejecta. The simulation is done 
with 600 cells at the lowest grid level, and we allow for 13 
levels on spatial grid —2.0 < x < 2.0 with a discontinuity 
at X = 0.0. In Fig |B.4| we plot the result at i = 0.4. In this 
test, we succeed to resolve the forward shock propagating 
to the right, the contact discontinuity, and the rarefaction 
wave propagating to the left (as the fluid in the left is hot). 
The left propagating rarefaction wave increases the normal 
speed to UxX ^ 0.879, and together with the transverse 
speed the Lorentz factor reaches the value 7 ^ 233. This 
high Lorentz factor in a thin shell between the contact dis- 
continuity and the tail of the rarefaction wave represents 
a very difficult numerical problem and requires the use of 
AMR. 



B.4. ID Riemann problem tests for Gamma Ray Bursts 

We now test the ability of the code to model the ultra- 
relativistic shocks between a relativistic cold flow with a 
cold external medium. This interaction induces a strong 
compression of the swept up matter from the external 
medium. We model three such cases, two with a Lorentz 
factor of 100 and a third with a Lorentz factor of 10. In all 
cases, we will find a reverse shock, contact discontinuity, 
and forward shock configuration. In these tests the den- 
sity of the fast flow is 10 times higher than the density of 
the stationary medium. All the parameters of the tests are 
given in Table [B4] These extreme conditions are represen- 
tative of GRB dynamics, and any code tailored for studying 
ultra-relativistic flows must be able to compute these situ- 
ations accurately. All these test are done with the hybrid 
version of HLLC. In all these tests there are all four regions 
that characterize the interaction between an outward mov- 
ing relativistic beam and the cold external medium: (1) the 
external medium at rest, (2) the shocked external medium, 
(3) part of the beam material which is shocked by the re- 
verse shock, (4) unshocked cold material of the beam. 
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Fig. B.5. GRB test case: a cold Lorentz factor 100 flow penetrates into cold static medium. The solid lines are the exact 
solution. Left: the density, center: Lorentz factor, right: effective polytropic index. 
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Fig. B.6. Second test case: a cold Lorentz factor 10 flow penetrates into cold static medium. The solid lines are the exact 
solution. Left: the density, center: Lorentz factor, right: effective polytropic index. 
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Fig. B.7. GRB test case with transverse velocity: a cold Lorentz factor 100 beam with also a small transverse velocity 
flow, penetrates into cold static medium. The solid lines are the exact solution. Left: the density, center: Lorentz factor, 
right: transverse velocity. 



In the first test, the simulation is done with 256 cells 
at the lowest grid level, and we allow for 12 levels on 
the spatial range —0.1 < a; < 10 with discontinuity at 
X = 0.0, the result at t = 9.5 is shown in Fig. |B.5[ The 
compression rate of the density by the forward shock is 
72 n2/ni_ = 72 (472 + 3) ~ 669.73 and the thermal energy 
reaches e2/'mp — (72 — 1) fT-2 ^ 616.86 ni as the Lorentz 



factor after the forward shock is 72 
where f = ^ 

77,1 



10 (Sari & Piran 19951. In the region 



V2 



12.57 



of the shocked swept-up matter, the polytropic index is 
4/3 as the thermal energy is higher than the mass energy. 
Moreover, the reverse shock is propagating with a relative 
Lorentz factor 73 ~ 2p/^ ^ 3.97, hence the matter is 

compressed according to n3/n4 ^ (473 + 3) ~ 18.9 (the 
exact obtained value is in fact somewhat lower) and the 
thermal energy enhanced es/rup = (73 — 1) 77.3 = 56.07 714. 
This implies that the state of the matter is also relativis- 
tic and the polytropic index r3 ~ 4/3. In this case both 
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Table B.l. GRB test parameters, the Lorentz factor and 
the transverse speed of left flow, and the density ratio be- 
tween left and right state. The fluid in the right is initially 
at rest. 



Test 


7L 


Pl/pr 




1 


100 


10 


0.0 


2 


10 


10 


0.0 


3 


100 


10 


0.01 



the forward shock and the reverse shock are relativistic. 
The numerical difliculty is to resolve the very thin region 
resulting from the strong compression of swept-up matter, 
which has an extent of /S.X2(t) = 0.0015t. It is easier to 
resolve the region between the reverse shock and the con- 
tact discontinuity, as the compression by the reverse shock 
is weaker, and this region extends over /S.X^{t) = 0.00315i. 

In a second test, a relativistic cold flow with a Lorentz 
factor 10 interacts with a cold external medium. This case 
corresponds more to the interaction of an AGN or micro- 
quasar jet with external medium than with a GRB case. 
The simulation is done with 256 cells at the lowest grid 
level, and we allow for 8 levels, the result dX t — 9.5 is 



shown in Fig. B.6 With the weak efficiency of the forward 
shock 72 ~ 3.79 to compress, the thermal energy of swept- 
up matter increases less and the effective polytropic index 
reaches a value Fefi = 1.356. The compression rate in this 
case is lower and hence the swept-up matter has an ex- 
tent of AX2{t) — 0.0163i, which allows to resolve easily 
the contact discontinuity and the forward shock. The re- 
verse shock is mildly relativistic with a relative Lorentz 
factor 73 ~ 1.53 such that the effective polytropic index of 
the shocked beam drops only to Fcff = 1.48. That makes 
the reverse shock less efficient to compress the beam matter 
(compared to the case of a shock with a lower polytropic in- 
dex F — 4/3), and the region between the reverse shock and 
contact discontinuity then extends over AX^{t) — 0.0225t. 

When we compare these two simulations, which differ 
in Lorentz factor of the beam, namely 100 versus 10, these 
tests show the dual numerical difficulty when the Lorentz 
factor increases: (1) the compression rate enhances and (2) 
the distance between the forward shock and contact discon- 
tinuity and between the reverse shock and contact discon- 
tinuity decreases. 

The last test is a variant of the first, with Lorentz fac- 
tor of the beam 100 and transverse velocity 774 j, = 0.01. 
Introducing a transverse velocity in the beam increases 
the numerical difficulty to resolve the contact discontinuity 
with a jump in both density and transverse speed. The sim- 
ulation is done with 256 cells at the lowest grid level, and 
we allow for 12 levels on the spatial range —0.1 < x < 10, 
the result at t — 9.5 is shown in Fig. |B.7[ As a result of 
the interaction, the transverse velocity in the beam mate- 
rial shocked by the reverse shock, increases to v^^^y — 0.015, 
which is a purely relativistic effect and depends both on 
the relativistic kinematics as on the state of the matter 
That mechanism could be of interest 
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in the refresii shocks for GRB afterglows, due to the ef- 
fect of the blastwave meeting a sudden (wind) ter mination 
shock in the Circ umBurst Medium density profile (Meliani' 
fc Keppens||2007| ). In more than 1 dimension, such mecha- 
nism will contribute to the shocked shell spreading with a 
speed higher than the sound speed. 



In these tests, it became clear that in the case when both 
the forward and reverse shock are relativistic, the constant 
polytropic index that should be used in simulations that do 
not have varying effective polytropic index is in fact 4/3. 
Even if the external medium in these tests is cold, the im- 
portant polytropic index is of the shocked fluids. In other 
cases when the density ratio 714/ni > 74 is very high, the re- 
verse shock is near-Newtonian and the effective polytropic 
index increases to approach 5/3 in the shocked medium be- 
tween reverse shock and contact discontinuity. Moreover, 
with a mildly-relativistic forward or mildly relativistic re- 
verse shock, as in the case with a beam Lorentz factor 
74 < 10, the effective polytropic index in swept up ISM 
matter is higher than 4/3 while in the shocked beam mat- 
ter the effective polytropic index can be mildly-relativistic 
as in Fig. |B.6| where Foff ~ 1.48. For all cases, it is more 
convenient to use the new EOS. 
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